/********************************************************************************
 *    Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH    *
 *                                                                              *
 *              This software is distributed under the terms of the             * 
 *              GNU Lesser General Public Licence (LGPL) version 3,             *  
 *                  copied verbatim in the file "LICENSE"                       *
 ********************************************************************************/
//*-- AUTHOR : Ilse Koenig
//*-- Modified : 11/11/2003 by Ilse Koenig
//*-- Modified : 16/05/99 by Ilse Koenig

/////////////////////////////////////////////////////////////
//
// FairGeoTrap
//
// class for the GEANT shape TRAP
//
// The technical coordinate system of a TRAP, which sits in
// CAVE and is not rotated, is the laboratory system.
// The y-axis points from the smaller side to the larger one.
// That's the same definitition as for a TRAP and different from
// the Geant or ROOT definition for a TRAP.
// Therefore a transformation is needed:
//              x-technical = - (x-Geant)
//              y-technical = - (y-Geant)
//              z-technical = z-Geant
// This is stored in the data element intrinsicRot which is
// created in the function calcVoluPosition(...)
//
/////////////////////////////////////////////////////////////
#include "FairGeoTrap.h"

#include "FairGeoTransform.h"           // for FairGeoTransform
#include "FairGeoVector.h"              // for FairGeoVector
#include "FairGeoVolume.h"              // for FairGeoVolume

#include <iosfwd>                       // for ostream
#include "TArrayD.h"                    // for TArrayD
#include "TMath.h"                      // for ATan, Pi, Sqrt, atan
#include "TMathBase.h"                  // for Abs
#include "TString.h"                    // for TString

#include <iostream>                     // for operator<<, cout, etc

using std::cout;

ClassImp(FairGeoTrap)

FairGeoTrap::FairGeoTrap()
  :FairGeoBasicShape(),
   intrinsicRot(FairGeoRotation())
{
  // constructor
  fName="TRAP";
  nPoints=8;
  nParam=11;
  param=new TArrayD(nParam);
  intrinsicRot.setElement(-1.,0);
  intrinsicRot.setElement(-1.,4);
}


FairGeoTrap::~FairGeoTrap()
{
  // destructor
  delete param;
  param=0;
  delete center;
  center=0;
  delete position;
  position=0;
}


TArrayD* FairGeoTrap::calcVoluParam(FairGeoVolume* volu)
{
  // calculates the parameters needed to create the shape
  const Double_t fac=20.;
  const Double_t raddeg=180./TMath::Pi();
  Double_t alpha, beta;
  FairGeoVector cb, ct, dc;
  for(Int_t i=0; i<4; i++) { cb+=*(volu->getPoint(i)); } // bottom plane
  for(Int_t j=4; j<8; j++) { ct+=*(volu->getPoint(j)); } // top plane
  dc=(ct-cb);
//  cout << dc <<" "<< ct <<" "<< cb <<" "<< endl;
  dc*=0.25;    // vector from bottom to top plane
  dc.setX(-dc(0));  // GEANT coordinate system
  dc.setY(-dc(1));
  param->AddAt(TMath::Abs(dc(2))/fac,0);
  alpha=TMath::ATan(TMath::Sqrt(dc(0)*dc(0)+dc(1)*dc(1))/dc(2))*raddeg;
  if (TMath::Abs(alpha)<0.0001) {
    alpha=0.0;
    beta=0.0;
  } else {
    if (TMath::Abs(dc(0))<0.0001) {
      if (dc(1)>0) { beta=90.0; }
      else { beta=270.0; }
    } else {
      beta=atan(dc(1)/dc(0))*raddeg;
      if (dc(0)<0) { beta=180.0 + beta; }
      if (beta<0) { beta=360.0 + beta; }
    }
  }
  param->AddAt(alpha,1);
  param->AddAt(beta,2);
  param->AddAt(((*(volu->getPoint(1)))(1)-(*(volu->getPoint(0)))(1))/fac,3);
  param->AddAt(((*(volu->getPoint(1)))(0)-(*(volu->getPoint(2)))(0))/fac,4);
  param->AddAt(((*(volu->getPoint(0)))(0)-(*(volu->getPoint(3)))(0))/fac,5);

  Double_t a=TMath::ATan(((*(volu->getPoint(1)))(0)
                          - (*(volu->getPoint(0)))(0)
                          + (*(volu->getPoint(2)))(0)
                          - (*(volu->getPoint(3)))(0))/40./param->At(3)
                        )*raddeg;
  if (TMath::Abs(a)<=0.0001) { param->AddAt(0.,6); }
  else { param->AddAt(a,6); }
  param->AddAt(((*(volu->getPoint(5)))(1)-(*(volu->getPoint(4)))(1))/fac,7);
  param->AddAt(((*(volu->getPoint(5)))(0)-(*(volu->getPoint(6)))(0))/fac,8);
  param->AddAt(((*(volu->getPoint(4)))(0)-(*(volu->getPoint(7)))(0))/fac,9);
  a=TMath::ATan(((*(volu->getPoint(5)))(0)
                 - (*(volu->getPoint(4)))(0)
                 + (*(volu->getPoint(6)))(0)
                 - (*(volu->getPoint(7)))(0))/40./param->At(7))*raddeg;
  if (TMath::Abs(a)<=0.0001) { param->AddAt(0.,10); }
  else { param->AddAt(a,10); }
// check if coplanar
  Double_t dx=(param->At(4) - param->At(5)) / param->At(3) * param->At(7) -
              (param->At(8) - param->At(9));
  if (TMath::Abs(dx)>=0.001) {
    cout << "top and bottom plane are not coplanar for shape TRAP\n";
    cout << "lenght in x-direction of top plane is changed\n";
    cout << "old values: " << param->At(8) << "  " << param->At(9) << "\n";
    param->AddAt(param->At(8) + dx/2.,8);
    param->AddAt(param->At(9) - dx/2.,9);
    cout << "new values: " << param->At(8) << "  " << param->At(9) << "\n";
  }
  return param;
}


void FairGeoTrap::calcVoluPosition(FairGeoVolume* volu,
                                   const FairGeoTransform& dTC,const FairGeoTransform& mTR)
{
  // calls the function posInMother(...) to calculate the position of the
  // volume in its mother
  Double_t t[3]= {0.,0.,0.};
  for(Int_t i=0; i<8; i++) { t[0]+=(*(volu->getPoint(i)))(0); }
  t[0]/=8.;
  t[1]=((*(volu->getPoint(1)))(1) + (*(volu->getPoint(0)))(1) +
        (*(volu->getPoint(5)))(1) + (*(volu->getPoint(4)))(1))/4.;
  t[2]=((*(volu->getPoint(4)))(2) + (*(volu->getPoint(0)))(2))/2.;
  center->setTransVector(t);
  center->setRotMatrix(intrinsicRot);
  posInMother(dTC,mTR);
}


